Robust ab initio calculation of condensed matter: transparent convergence through 

semicardinal multiresolution analysis 
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We present the first wavelet-based all-electron density-functional calculations to include gradi- 
ent corrections and the first in a solid. Direct comparison shows this approach to be unique in 
providing systematic "transparent" convergence, convergence with a priori prediction of errors, to 
beyond chemical (millihartree) accuracy. The method is ideal for exploration of materials under 
novel conditions where there is little experience with how traditional methods perform and for the 
development and use of chemically accurate density functionals, which demand reliable access to 
such precision. 

PACS numbers: 71.15.Ap, 71.15.Mb, 71.15.Nc 
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Over the last several decades, the ah initio density- 
functional approach, which replaces the many-body wave 
function with a set of single-particle "Kohn-Sham" or- 
bitals moving in an effective potential^], has opened to 
first-principles study a diverse array of condensed mat- 
ter phenomena ranging from plasticity, diffusion and sur- 
face reconstruction to melting and chemical reactions 
While density-functional theory is exact in principle, the 
true form of the effective potential as a functional of 
the orbitals is unknown. Current research is pushing 
available approximations to this functional toward chem- 
ical accuracy 3], typically defined as predicting bond- 
breaking energies to 1 kcal/mol — 1.6 millihartree. The 
prime motivation for this precision is the ability to pre- 
dict rates of microscopic processes at room temperature, 
where the relevant energy scale is k^T k, 1 millihartree. 

Unfortunately, current feasible representations for the 
Kohn-Sham orbitals, such as the plane-wave pseudopo- 
tential approach^ or the atomic sphere methods'^!, are 
not easily improved systematically to millihartree accu- 
racy. For instance, although plane-wave bases converge 
systematically with basis size, pseudopotential calcula- 
tions require significant experience in the construction 
of reliable potentials 0]- And, although atomic sphere 
methods use the Coulomb potential of the nuclei directly, 
they remain also an art involving many parameters and 
requiring significant expertise to ensure millihartree ac- 
curacy consistently (5(1 . As testimony to these difficulties, 
it is not uncommon to find in the literature disagreements 
over the predictions of density-functional theory. 

To illustrate this. Table HI presents the lattice con- 
stant, cohesive energy and bulk modulus of a very simple 
solid, MgO in its rock-salt structure, both as measured in 
experiments!^ ^"■c^ and as calculated within the lo- 
cal spin-density approximation (LSDA) when using yari- 
ous representations for the Kohn-Sham orbitals |9l llOl [ill 
[12 . 13] . No two methods in the table agree consistently 
on the predictions of LSDA. For each property, the spread 
in predictions is on the order of the discrepancy from the 





a [A] 


Ecoh [eV] 


B [Mbar] 


Experiment f6, 7, 8] 


4.21 


10.33 


1.55 - 1.62 


LMTOfg] 


4.09 


10.67 


1.71 


FP-LMTOflO] 


4.16 






FP-LAPWfllj 


4.167 




1.72 


pseudopotentialfl2] 


4.191 


9.96(+0.5) 


1.46 


pseudopotentialflS] 


4.125 


11.80 


1.56 



TABLE I: Structural properties of MgO 



experiment. Without direct estimation of the errors, it 
is difficult to judge how much of each discrepancy is in- 
herent in LSDA or due to biases built into the different 
representations. Hence, the use of existing methods is 
questionable when exploring materials under novel con- 
ditions where there is little experience with how such 
methods perform, as in the study of materials at geolog- 
ical pressures under which atomic cores begin to overlap. 
(There is recent interest in ah initio study of MgO under 
such conditions!!^.) Finally, lack of consistent access to 
millihartree precision hampers development and eventual 
use of chemically accurate functionals. 

The purpose of this letter is to demonstrate for the 
first time that a new general representation, a wavelet- 
like multiresolution analysis, provides for solid-state elec- 
tronic structure calculations an unprecedented level of 
precision and a new capability for transparent conver- 
gence, systematic convergence with an extremely simple 
and predictable scaling for the errors. We also demon- 
strate the first wavelet calculations to employ generalized 
gradient approximations (GGAs). We find GGAs to fit 
seamlessly within our approach, without special concerns 
such as the discontinuities at the sphere boundaries which 
arise in the atomic sphere methods|l5j. 

To date, the application of multiresolution analysis 
to ah initio calculations has been limited solely to very 
simple systems, such as sin gle- electron atoms [l^. the 
diatomic hydrogen molecule [l7f. 



diatomic oxygen us- 
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ing pseudopotentials[l8l| . purely electrostatic problems 
without electronic structure or all-electron calcu- 
lations of atoms li^l- Only recently have all-electron 
wavelet calculations of small molecules appeared in the 
literature^21j 2^. Here, new techniques 23, 24j enable 
us to present the first such calculations in solids and the 
first to include gradient corrected density functionals. 

Multiresolution analysis — A multiresolution analy- 
sis consists of a basis set of spatially localized func- 
tions which describe fluctuations on a hierarchy of length 
scales, each separated by a factor of two, with the basis 
functions representing each of these levels of resolution 
organized on simple rectilinear grids. (Ref. provides 
a detailed review.) The central, nontrivial mathematical 
result of multiresolution analysis is that, with appropri- 
ately chosen basis functions, this multilevel description 
is mathematically equivalent to a uniforrn grid of basis 
functions on the finest level of resolution 25, '^E]. This 
key result allows for a priori knowledge of convergence. 
In particular, the use of third-order interpolating basis 
functions in the present work implies that all errors scale 
as the fourth power of the spacing on the finest level, a 
factor of sixteen for each additional level of resolution. 

The superiority of the multiresolution analysis over the 
uniform representation comes from the unique way in 
which the analysis represents information. Because fine- 
scale coeHicients carry information only about high spa- 
tial frequencies, they drop rapidly to zero with distance 
away from the nuclei |2lll24l|. Thus, restriction of the ba- 
sis, elimination of functions from finer levels far from the 
nuclei, has controllably negligible effect on the outcome 
of the calculation. In practice, sufficient functions can be 
restricted from the basis so that the basis size and work- 
load grow linearly with the number of levels of resolution, 
making convergence exponential with basis size |20l[2l l|. 

The primary reason for the limitation of wavelet cal- 
culations to very simple systems in the past has been 
the lack of efficient algorithms suited to the solution of 
non-linear partial differential equations. We recently in- 
troduced new algorithms 0,01 which are faster than the 
approaches used in the older wavelet electronic-structure 
works by some three to four orders of magnitude. With 
these methods, all-electron calculations now require an 
effort of the same order of magnitude as their pseudopo- 
tential counterpartsjl^, opening the possibility of all- 
electron wavelet calculations of non-trivial systems. 

Given a multiresolution representation, our density- 
functional calculations proceed by strai ghtf orward ex- 
pansion of the Kohn-Sham Lagrangian'24f,'27] or, equiva- 
lently, energy functional Jj, ernploying nuclear potentials 
constructed as described in [22| and employing either 
the Vosko, Wilk and Nussair parameterization (VWN) 
of LSDA 28, 29J or the Perdew, Burke and Ernzerhof pa- 
rameterization (PBE96) of GGAH^]. The calculations 
then locate the stationary point of the functional us- 
ing standard preconditioned conjugate gradient methods 
with analytic continuation as in j3ll | to maintain the or- 
thonormality constraints among the Kohn-Sham orbitals. 



Refs. mm give the full details of our implementa- 
tion, with the exception of four extensions which proved 
critical in carrying out calculations in solids to high pre- 
cision: (1) treatment of boundary conditions for Bloch 
states, (2) symmetrization of the electron density, (3) ex- 
pansion of the electron density on higher resolution grids, 
and (4) the extension of our approach to include gradient 
corrections to the exchange-correlation energy. 

(1) In a finite multiresolution basis, the common prac- 
tice of expanding the periodic parts Unk {t) , rather than 
the full Bloch orbitals ipnkir) = e*'^''u„fc(r), does not 
ensure extensivity, that decreasing cell size while corre- 
spondingly increasing Brillouin-zone sampling results in 
identical total energies. However, because wavelet bases 
share the translational symmetries of the lattice, their 
expansion coefficients satisfy a discrete Bloch's theorem, 
Cx+R = e^'^^Cx, where Cx and R are, respectively, the ex- 
pansion coefficient for the basis function centered at x 
and any lattice vector. We thus implement Bloch states 
by storing the coefficients Cx for each orbital tpnk in a 
single representative cell and producing all needed coeffi- 
cients Cx+R by multiplying the corresponding coefficient 
Cx by e*'^^ ( "twisted boundary conditions" ) . 

(2) For k-point sampling, we use the scheme of 
Monkhorst and Pack jsj with symmetry folding of the 
Brillioun zone. For this folding to be exact, the sym- 
metrized electron density must be expandable in a mul- 
tiresolution basis respecting the symmetries of the crys- 
tal. Accordingly, we symmetrize the density with the 
operator s = IPJSIPJ, composed from the usual 
real-space symmetrization operator 5, a projector onto 
basis functions which have all symmetry partners P, 
and the forward I and inverse J wavelet transforms of 
[2^. Moreover, differentiation of the total energy shows 
that the self-consistent potential requires a further sym- 
metrization with the operator s\ a different operator 
from s because J ^l'^ for non-orthogonal bases. 

(3) In our older approach, which samples real-space 
quantities only at the centers of basis functions |22j|. errors 
in approximate evaluation of the total energy dominate 
errors from incomplete representation of the orbitals. Ac- 
cordingly, we now save considerable computational effort 
by expanding the orbitals in a basis of exactly half the 
spatial resolution of that representing the other quanti- 
ties. This proves critical to the present calculations and 
requires the introduction of new transform operators to 
be described in depth in a forthcoming publication. 

(4) Gradient corrections implement simply and natu- 
rally within our framework. Using the notation estab- 
lished in [2^ , the exchange-correlation energy 

Exc = j d^x n{x)cxc {dxin{x),dx2n{x), dx3n{x) , n{x)) 

becomes 

Exc - (Jn)^ OJexcCDiJn, V2Jn, VsJn, n), (1) 

where n is a vector of values of the electron density at 
sampling points at the centers of the basis functions, the 
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Vi are matrices of values of the derivative of each 
basis function at each such sample point, and O is the 
matrix of overlaps among basis functions. The exchange- 
correlation potential then follows directly by differentiat- 
ing 1^ with respect to n, which introduces no new oper- 
ators other than the Hermitian conjugates of the P;. 

Transparent convergence — To illustrate the feasibil- 
ity of larger calculations in complex systems, we study 
an eight atom cubic supercell of MgO. Without pseudiz- 
ing or treating differently the core states or core regions, 
our results converge directly to the true LSDA/GGA pre- 
dictions as a function of only three parameters: number 
of iterations in the self-consistent solution of the Kohn- 
Sham equations, the number of k-points in the sampling 
of the Brillouin zone, and the size of the multiresolution 
basis. The remaining paragraphs of this section address 
these three parameters one by one. 

Despite the broad range of length scales in the calcu- 
lation, the convergence of conjugate-gradient minimiza- 
tion is extremely good (0.12 digits/iteration), even when 
compared to that of plane-wave pseudopotential calcu- 
lations. As discuss, this results from the use of 
our specific non-orthogonal basis with a simple diago- 
nal preconditioner. This approach achieves millihartree 
precision (corresponding to six significant figures) within 
just forty iterations of starting from randomized atomic 
wave functions. We typically ran eighty iterations. 

Exploring convergence with respect to Brillioun zone 
sampling (with a multiresolution basis somewhat smaller 
than that which the calculations below employ) we find 
that that eight k-points, which reduce to one special point 
at [0.25 0.25 0.25] in the irreducible wedge, suffice to con- 
verge the energy to 0.3 millihartree per chemical unit, 
better than chemical accuracy. 

A key, novel result of this work is the simple form of the 
convergence to the full all-electron result with increasing 
basis set and the high precisions which this allows us to 
reach. Figuren]shows for both LSDA and GGA the total 
energy per MgO chemical unit as a function of the fourth 
power of the spacing on the finest level for the multireso- 
lution grids in Table ITU when truncated at four, five and 
six levels of refinement. The results demonstrate that the 
basis convergence of GGA is almost identical that that 
of LDA, with no adverse effects from the presence of sec- 
ond derivatives of the density in the exchange-correlation 
potential. The data clearly exhibit the a priori expected 
quartic convergence, thus demonstrating that the restric- 
tion of the basis has negligible effect and that the cal- 
culation has entered the asymptotic convergence regime 
where there are no hidden convergence "shoulders" . The 
quality of the linear fit empowers extrapolation to infinite 
resolution (stars in the figure insets) with an error of only 
~14 TTizcrohartree, far below the tolerances described for 
the most accurate of the standard approaches^. We 
also see that, even without extrapolation, six levels of 
refinement suffice to give the energy to within 0.5 milli- 
hartree per chemical unit. Thus, the present approach 
gives simple, transparent knowledge of the precision of 
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FIG. 1: Convergence of total energy per MgO chemical unit 
with basis size: LSDA (left panel), GGA (right panel). Insets 
are comparable in size to square symbols in main plots. 



Level 







44, 44, 44 




Mg 





1 


24, 24, 24 


24, 24, 24 


2 


24, 24, 24 


24, 24, 24 


3 


24, 24, 24 


28, 28, 28 


4 


28, 28, 28 


28, 28, 28 


5 


32, 32, 32 


32, 32, 32 



TABLE II: Restriction employed in the present work: dimen- 
sion {N:^, Ny, Nz) of the cubic grids of basis functions on the 
top coarse periodic level (0), and lower levels (1-5) centered 
on magnesium and oxygen nuclei. 



the calculation at any stage. We believe the simplicity 
of the convergence of multiresolution analysis to the all- 
electron result in a practical calculation to be unique in 
the field of electronic structure. 

Results — Table IIIII summarizes the predictions of 
the LSDA-VWN and GGA-PBE96 parameterizations of 
density-functional theory, when computed with one spe- 
cial k-point and the six-level restriction of the multireso- 
lution analysis from Table Ull which the discussion above 
establishes to give millihartree precision. Note that, to 
allow direct comparison with the values available in the 
literature, the table lists the value of the binding en- 
ergy at the experimental lattice constant. Comparison 
of Tabled with Table imi shows that, of the traditional 
methods, only FP-LAPW consistently reproduces the 
fully converged LSDA results to chemical accuracy. Fur- 
ther, it is now possible to make a definitive determina- 
tion of the relative transferability among pseudopoten- 
tials. Finally, we can now judge unambiguously the pre- 
dictions of LSDA and GGA for MgO. As generally the 
case for LSDA, the lattice constant is within a few per- 
cent (-1.5%) of the experiment, the system is slightly 
over bound (-1-14%), and the bulk modulus is in error 
by several percent (5-10%, depending upon to which ex- 
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a [A] 


Ecoh [oV] 


B [Mbar] 


LSDA 


4.161±0.003 


11.80±0.03 


1.71±0.03 


GGA 


4.221±0.004 


11.90±0.03 


1.64±0.03 






10.4+ 





+ Including estimate for spin-polarization energy of atomic 
oxygen based on LSDA/LDA difference of 1.5 eV. 

TABLE 111; Multiresolution analysis calculation of structural 
properties of MgO within LSDA and GGA. 

periment we compare). For GGA relative to LSDA, we 
find the typical expansion of the lattice constant, cor- 
rection of the over-binding (when including atomic spin- 
polarization energies) and decrease in the bulk modulus. 
For MgO, GGA improves these results dramatically. 

Conclusions — We report the first solid-state all- 
electron calculations within a wavelet-like multiresolu- 
tion analysis and the first multiresolution calculations to 
include gradient corrections. The results give the first 
independent verification that, of the traditional meth- 
ods, FP-LAPW is by far the most accurate and gives 
the first unambiguous demonstration that GGA outper- 



forms LSDA for MgO. We also demonstrate multireso- 
lution analysis to be unique among electronic-structure 
representations in exhibiting extremely simple and pre- 
dictable convergence to far beyond millihartree accuracy. 
This approach, therefore, is ideal for novel problems such 
as the study of the impact of core polarization effects on 
electron energy loss and X-ray absorption spectra and 
the investigation of materials under geological pressures, 
applications where the use of the traditional, less system- 
atic approaches remains more an art than a science. 
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